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Abstract 

We show that the recent hierarchy of semidefinite programming relaxations based 

on non-commutative polynomial optimization and reduced density matrix variational 

o ' methods exhibits an interesting paradox when applied to the bosonic case: even though 

it can be rigorously proven that the hierarchy collapses after the first step, numerical 

C . implementations of higher order steps generate a sequence of improving lower bounds 

^ , that converges to the optimal solution. We analyze this effect and compare it with 

O^l similar behavior observed in implementations of semidefinite programming relaxations 

for commutative polynomial minimization. We conclude that the method converges due 

CN . to the rounding errors occurring during the execution of the numerical program, and 

show that convergence is lost as soon as computer precision is incremented. We support 

[■*-. ■ this conclusion by proving that for any element p of a Weyl algebra which is non- 

t^ , negative in the Schrodinger representation there exists another element p arbitrarily 

close to p that admits a sum of squares decomposition. 
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1 Introduction 

Computing the energy spectrum of a finite set of indistinguishable particles subject to a 
given potential is a standard problem appearing in many branches of physics, e.g., in quan- 
tum chemistry, atomic physics, or condensed matter physics. Although traditionally the 
main approaches to this problem have been variational [T], in the last decade it has been 
attacked with success by means of semidefinite programming (SDP) formulations [21 [3l H] 
of the constraints on second-order reduced density matrices proposed in [5l El [7] (the so- 
called 2-RDM method). These SDP methods can be viewed as particular instances of a 
more general non- commutative polynom,ial optimization approach [Si [9], which extends to 
the non-commutative setting the method developed by Lasserre [TD] and Parrilo pLlj for 
scalar polynomial optimization. Roughly speaking, a non- commutative optimization prob- 
lem consists in finding the minimal eigenvalue of a hermitian polynomial of non-commutative 
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operators. To solve such problems, one can define a hierarchy of SDP relaxations, each of 
which corresponds to finding a sum of squares decomposition of the polynomial to be mini- 
mized which provides a lower bound p^ on the optimal solution p* of the original problem, 
with p^ < p^ < ... < p*. This approach reduces to the 2-RDM method when applied to 
fermionic systems, but more generally is also highly successful, e.g., to characterize the set 
of quantum correlations in quantum information science [12]. Recently, modifications of this 
algorithm exploiting translational invariance have been proposed independently by Hiibener 
and Barthel [13] and Baumgratz and Plenio [H] as an alternative to variational techniques 
in condensed matter physics. In [151 [131 E]) it was suggested to apply such SDP methods 
to compute the ground-state energy of bosonic systems, i.e., to find the minimal eigenvalue 
of Weyl polynomials. 

In this article, we point out that any computer implementation of the SDP hierarchies [H 
[9] to bosonic systems will exhibit the non-commutative analog of an effect already observed 
in similar algorithms for commutative polynomial minimization [16], [17]. On one hand, it 
can be proven that any relaxation beyond the first one will not provide better lower bounds 
on p*, i.e., p^ = p^ for all k > 1. On the other hand, numerically it is observed that the 
bounds p^,p^, ... output by the computer form an increasing sequence, with limfc_>ooP^ = P*- 

We will show that the resolution of this "mathematical paradox" follows the same lines as 
the commutative one. Namely, even though there exist positive Weyl polynomials p that do 
not admit a sum of squares decomposition, for any such polynomial there exists an arbitrarily 
small perturbation p ^ p = p + eg such that p can be decomposed as a sum of squares of 
Weyl polynomials. The rounding errors introduced by the computer while executing the 
algorithm correspond to such a perturbation, and so numerical implementations of the SDP 
method converges to the correct answer of the problem. This is, therefore, an example of a 
numerical method that converges, not in spite of rounding errors, but because of them. 

This article is structured as follows. In Section [21 we provide the basic definitions and 
facts about Weyl algebras that are used in the remaining of the text. In Section [3] we 
will describe the SDP method and illustrate the paradox with a numerical example. The 
resolution of the paradox will be given in Section IU where we will prove that any positive 
Weyl polynomial can be perturbed to a sum of squares of polynomials. Finally, in Section [5] 
we will present our conclusions. 

2 Definitions and basic results on Weyl algebras 

2.1 Weyl algebras and Weyl polynomials 

A Weyl algebra Wn is a *- algebra with 2n generators a = (ai, 02, ..., a„) and a* = (a^, 03, ..., a: 
satisfying the canonical commutation relations (OCRs): 



' ""n/ 



[ai,aj] = [a*, a*] = 0, [ai,a*] = 6ij, for all i,j G {l,...,n}. (1) 

An element p of VV„ is thus a linear combination (with complex coefficients) of words in the 
2n letters a and a*. The words and elements of W„ can also be viewed as monomials and 
polynomials, respectively, in the 2n variables a and a*. 



Using the CCRs, any element p of a Weyl algebra can be brought to the normal form 

p = 5^p,-Ka*)V, (2) 

s,t 

where s = (si, . . . , Sn),i = (ti, . . . , t„) E N"", a* denotes the monomial a* = YYi=i ^^i ^^"^ 
similarly {a*y = YYi=i^V'j ^^^ where Ps^t £ C are complex coefficients. For instance, the 
elements aia* and 0102030^, expressed in normal form, are 1 + a^ai and 1 + a^ai + 0203 + 
0^030102, respectively. We will later show that the decomposition ([2]) is unique for each 
p G W„. To check whether two different polynomials p,p' in the variables a, a* represent the 
same element p = p' E Wn it is therefore enough to write p, p' in normal form. 

This last observation suggests a natural norm in W„: let p G >V„, and let ([2]) be its 
normal decomposition. Then, we define the /i-norm of p as 

hip)^Yl\P^rt\- (3) 

s,t 

It is also useful to distinguish the elements of >V„ by resorting to the concept of degree. We 
say that the degree of an element p of W„ is deg(p) = max{||s||i + ||t||i : Psi 7^ 0} for p 
expressed in normal form ([2]). It is easy to see that the degree of any monomial s = Y[k=i ^k 
with Sk G {aj,aj}"=i is equal to deg(s) = \s\ = d. 

Given a monomial s, we denote by |s| its anti-normal ordering, that is, the monomial s' 
that results when we reorder the letters appearing in the expression of s in such a way that 
all the letters ai, 03, ..., a„ end up on the left and all the letters (a*, Og, ..., a*) end up on the 
right. For example, Xa*2{a\YaiX = $0*010^02$ = Oi (0^)^03. 

In this article, we will be mainly concerned with hermitian elements of VV„, i.e., those 
polynomials p such that p = p*. If p is decomposed as in ([2]), the hermiticity condition thus 
translates as Ps^t = p|g, where p|g denotes the complex conjugate of the number pj g. 

2.2 The Schrodinger representation 

Let TT be a mapping n : >V„ — t- LiTi), for some Hilbert space "H, where LiTi) denotes the 
space of linear (not necessarily bounded) operators of H. We say that vr is a ^-representation 
of Wn if and only if 

1. 7r(l) = I^. 

2. TT{pq) = -n:{p)-n:{q), 7r(p + q) = 'n{p) + Tx{q) for all p, g G >V„. 

3. ■k{p)* = 7r(p*), for all p G W„. 

We now show that W„ admits a representation. For this, let H he a separable Hilbert space, 
and let {|s), s G N} be an orthonormal basis for H, which we call the number basis. If we 
denote by a G L{H) the linear operator defined by 

., , _ f for s = 0, 

\ y/s\s — 1) otherwise. 



then its adjoint a* satisfies 



~a*\s) = V7+l\s + l), (4) 

and so it can be verified that 

[a,a*]|s) = |s), for all s G N. (5) 

Defining l-i = iJ®"", we can then build a representation tt^ : VV„ — > -^("H) for the Weyl algebra 
W„ through 

7r5(afc)=I®'=-i®S®r"-^ (6) 

This representation of VV„ is known as the Schrodinger (or Fock) representation. From now 
on, we always refer to this representation and write n for ns for simplicity. 

2.3 Weyl polynomial minimization 

The Schrodinger representation admits a clear physical interpretation: given a set of n one- 
dimensional particles, it associates to each particle k G {l,...,n} a pair of creation and 
annihilation operators 7r(afc), 7r(a^). The operators describing the position and momentum 
of particle k along the real line are then given, respectively, by n^Xk) = 7r( °''^'' ) and n^pk) = 

T^i-^h^)- If these particles are subject to a potential y(a;i, ...,a:„), the energy operator of 
the system, in non-relativistic approximation, will be given by 

<E) = ^\y^^^+V{x,,...,xS\. (7) 

where mi G M"*" is the mass of particle i. In particular, the minimum energy of the system 
will be given by 

K,,{E) = mimniEM : |0) G S, (0|0) = 1}, (8) 

where S is the Schwartz space, that is, the set of states |0) (i.e., vectors of H) satisfying 
IK(p)0)ll < '^ fo'^ ^11 P ^ ^n- Note that the minimization over {(j)\E\(j)) makes sense, 
because the energy is an hermitian operator, E = E*, and consequently, {(f)\E\(f)) G M is a 
real quantity for all |0) E S. 

The case where E is a polynomial in the variable Xi,pj - or equivalently in the vari- 
ables Oj, a* - is particularly important (it includes for instance the case where the potential 
V{xi, ...,Xn) is Taylor expanded around some equilibrium position). This motivates the 
following generic Weyl polynomial optimization problem 

Ainf(p) = inf {(0|7r(p)|0) : |0) G S, (0|0) = 1}, (9) 

for an arbitrary hermitian polynomial p G VV„. Note that, alternatively, we can write 

Ai„f(p) = sup {A G M : 7r(p) - A > 0}, (10) 

where positivity is understood in the Schwartz space. This reformulation of the problem will 
be used in the next section. 



2.4 Coherent states 

An interesting subset of % is constituted by the coherent states. For any a G C, denote by 
\a) E H the normahzed state 

ap ■^^ — ^ ot^ 



seF * ^■ 



Then, a coherent state in l-i is any state of the form \a) = ®"=i|ai), for any a G C". Coherent 
states are important because they are simultaneous eigenstates of the annihilation operators 
{7r(afc) : k = 1, ..., n}. This follows from the easily verified identity d\a) = a\a) for all a G C. 
As a result, for any normally-ordered polynomial p G Wn, we have that 

(Q;|7r(p)|a) = y^Ps,tOi*^a^ = p{a,a*). (12) 

As an application of the coherent states, let us show that the decomposition ([2]) is unique, 
or, equivalently, that admits a unique representation. Suppose thus that 

= $^Paa*V, (13) 

s,t 

for some coefficients Ps^t G C. Using the Schrodinger representation, we have that 

= (a|7r(0)|a) = ^pg,i-a*V for all a G C". (14) 

s,t 

The right-hand side is a polynomial in the complex variables a, a*, and it can only be equal 
to zero for all values of a if ps^t = for all s, i. 

3 SOS decompositions, the SDP hierarchy, and the 
paradox 

Given a hermitian polynomial p G W„, we say that p admits a sum-of-squares (SOS) decom- 
position if there exist polynomials fi G >V„ such that 



p 



J2f*f^- (15) 



We denote E^ the set of all such polynomials. It is clear that if p G S^, then 7r(p) > 0, since, 
for any |0) G 5, 

(0i7r(p)i0) = Y.{<pw:Mm<p) > 0. (16) 

i 

However, the opposite implication is not true, not even in Wi. Indeed, as noted by Schmiidgen 
[18j, the family of polynomials Pe = {alai — l)(aia^ — 2) + e satisfies vr(pj > 0, for e > 0, 
but nevertheless p^ ^ S^, for e < i. 



Given a p G VV„, a possible scheme for finding a lower bound A* on the solution Ainf(p) 
of firUjl is thus to solve the problem 

A* = sup{AGM:p-AG S^}. (17) 

This principle is the one behind the polynomial minimization algorithms developed by 
Lasserre [10] and Parrilo [H], and their non-commutative analogue [HI [9]. Such algorithms 
work by searching for SOS decompositions of p — A with some degree constraint. Applied to 
the Weyl minimization problem, this results in the following sequence of programs: 

A'= = sup{AgM:j9-AgS^}. (18) 

Here k is an integer such that 2k > deg(p) and indexing the successive programs in the 
sequence and S^ is the set of polynomials which admit a decomposition of the form flT^ 
with deg(/j) < k, for all i. Each of these problems is a semidefinite program, as one can 
check that 

A'^ = max{A eR:p-\ = {w^yZw\ Z > 0}, (19) 

where w'' is a vector whose components are the normally ordered monomials of degree < k. 
Clearly, A'^ < A'^"'"^ and liniTv^oo A'^ = A*. The programs flTSl) . (1191) thus form a converging 
hierarchy of SDP relaxations for the problem (I17p . Supplemented with a boundedness con- 
dition (that is not satisfied in the present case of Weyl polynomials), it can be further be 
shown that this hierarchy necessarily converges to the optimal solution of the problem fITOl) 
[HIS], as problems (ITTj) and fITOl) then turn out to be equivalent [19]. 

Let L denote an arbitrary functional on the Weyl algebra, i.e., L : VV„ — ?■ C Then, the 
dual of problems f lTS]) . flT^ can be shown to be 

/ = min{L(p) : L(l) = 1, L{qq*) > for all q G >V„, deg(g) < k}, (20) 

or explicitly in SDP form 

/ = mm{J2p-s,ty-s,t : 2/0,0 = 1, M^iy) > 0}, (21) 

s,t 

In this last formulation, Pst are the coefficients of p in normal form ([2]), y = {^jf : ||s + t||i < 
2k} C C are the optimization variableqj, and Mk{y) is the moment matrix of order k, a 
matrix whose rows and columns are indexed by pairs of vectors (s, t) G N" x N*^, with 
11^ + t||i < k, and with entries defined by 

[^fc(?/)](.,-),(.,.) = E^'?.^-^5.^-' (22) 

q,r 

where {cq^f} are the normal form coefficients of the monomial {a*'^a^)*a*^a^, i.e., 

{a*~'a'ya*''a^ = ^ c^-^.-a* V. (23) 

q,f 

The next lemma shows that problems (TTHj) . (TT9l) and (120|) . (12T1) are, in fact, equivalent: 



'^Note that, ii pg.i G R,Vs,t, one can assume {j/j,* ■ \\s + t\\i < 2k} C 



Lemma 1. // there exists a feasible point of program / fTPj) . then ii^ = \^ . 

Proof. By Sylvester's criterion |2Qj, it just suffices to show that problem f l20|) admits a strictly 
feasible point, i.e., that there exists a functional L such that -^(1) = 1 and L{pp*) > 0, for 
all p y^ 0. Now, consider the functional L{p) = tr(n7r(p)), with 

n = j^Jdae-\-\'\a){a\. (24) 

This functional satisfies L{1) = 1. Also, for any non-zero polynomial p ^ E Wn we have 
that 

L{p*p) = ^y"dae-l^l^«|7r(p»|«)>^|dae-l^l^|(«|vrb)|«)r 

^ ^dae-l"l'|p(a,a*)|2>0. (25) 



(27r)' 

The last inequality comes from the fact that p{a, a') ^ and that the integration takes place 
in all C". D 

Even though the above SDP hierarchy cannot be guaranteed to converge to the optimal 
value Ainf(p) of ( ITU]) , every SDP step provides a lower-bound on Ainf(p). Based on the 
successful applications of this SDP hierarchy to fermionic systems and quantum correlations, 
where in practice very good lower bounds are obtained after only a few SDP relaxations, one 
could expect a good overall performance also in the context of Weyl polynomials. However, 
as the next result shows, no improvement over the first lower-bound can be obtained by 
considering higher steps in the hierarchy. 

Lemma 2. Let p eT? . Then, p E T^l for 2k < deg{p). 

Proof. By hypothesis, p G S^, and thus there exist polynomials /, such that (^ holds. We 
will show that all such polynomials satisfy 2deg(/j) < deg(p). 
Suppose, on the contrary, that 

with 2deg(/) > deg(p). Then, 

p{a,a*) = {a\7r{p)\a) > {a\7r{f*)7T{f)\a) > \{a\iT{f)\a)\' = \f{a,a*)\'. (27) 

Now, denote by LT{f) the leading terms of /, i.e., 

LT{f){a,a*)= Yl /^-*V)*«', (28) 

||s+f!|i=deg(/) 

and choose /3 G C" such that LT{f){f3,f3*) = c ^ 0. Then it is straightforward that 
|/(r/3,r/3*)| = ©(r^^S^^)). However, \p{r^,r^*)\ < 0(r^^g(p)). For (ET]) to hold, we must 
thus have that 2deg(/) < deg(p). D 
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Figure 1: Bounds on Ainf(-Em) as a function of m. The upper curve represents an upper 
bound obtained through variational techniques. The lower curve represents the lower bound 
A^ corresponding to the first relaxation of the problem. 

What Lemma [2] shows is that, for any polynomial p, the sequence of values A'^'', A'^"^"^^, ..., 
with fco = [deg(p)/2] is constant and equal to A*. In other words: the first SDP relaxation 
of the problem already provides the best approximation to Ainf(p) attainable with SOS 
decompositions. 

How does such an approximation perform? Consider the uniparametric family of one- 
dimensional hamiltonians {Em '■ ^^ ^ I^}) with 



Er, 
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-p + vnx + X , 



(29) 



with x = (a + a*)/ \/2^ p = (a — a*)/i\/2. Note that, for m < 0, Em corresponds to the 
interesting double-well potential. Figure [T] shows a plot of A^ as a function of m, together 
with an upper bound on Em obtained through variational methods. We used the solver 
SDPT3-4.0 [2T] and the MATLAB package YALMIP [22] to carry out the SDP calculations. 
It is clear that, as soon as m < 0, the approximation given by A^ becomes worse. 

From our discussions above, it follows that A^ = A^ = A^ = ..., and so A^ represents the 
best lower bound on Ainf(em) achievable with the SDP hierarchy. Figure [2] shows, however, 
that such is not the case. Indeed, we see that subsequent relaxations of the problem return 
lower bounds which are closer and closer to the variational upper bound, until, at A^, both 
bounds become practically indistinguishable. What is happening? 
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Figure 2: Lower-bounds on Ami (-Em) as a function of vn corresponding to the second, third, 
fourth and fifth SDP relaxations, in ascending order. The black line represents the variational 
estimate, quite close to A^ (red line). 



4 Resolution of the paradox 

As mentioned in the introduction, the above paradox is not new in commutative polynomial 
optimization: indeed, Henrion and Lasserre [16] noticed that the numerical implementation 
of their SDP algorithm for polynomial minimization returned the optimal value of the 2- 
dimensional Motzkin polynomial, instead of its SOS value (— oo). Lasserre successfully solved 
this paradox by proving that any commutative positive polynomial can be approximated ar- 
bitrarily well by an SOS decomposition [23j- The accepted resolution of the paradox was that 
the rounding errors occurring during the numerical computations perturbed the polynomial 
p to be minimized to another one of higher degree admitting an SOS decomposition |17] . 

In this Section we will prove a non-commutative analog of this result, namely, that 
Weyl polynomials which are positive semidefinite in the Schrodinger representation can be 
perturbed to a higher degree polynomial in Y? . This is formally stated in the following 
Theorem: 

Theorem 3. Let "p he an element ofWn such that tt{p) > 0. Then, for any e,S > 0, there 
exists a polynomial p ^T? such that li{p — p) < e and \\m{{p) — Ainf(p)| < S. 

The proof of this theorem follows straightforwardly from the next three lemmas. In these 
lemmas, the constant c is arbitrary but fixed to be c > 2. 



Lemma 4. Let p G VV„, be a polynomial such that 7i{p) > 0. Let 

(n~l)\ 

\\t\\i<r 



^^ A^ c\\4i(n+ ||t||i - 1)! ^ ^ ' ^ ^ 



where the sum runs over all vectors t ^N^ of length \\t\\i less or equal than r. Then for any 
e > 0, there exists some tq such that for all r > tq, 

f=p + eg:eT.\ (31) 

Lemma 5. Let p^ he defined as in Lemma^ Then 

h{p~f)<e-^, (32) 

for any r > 0. 

Lemma 6. Let p be a hermitian polynomial in Wn such that Xm{{p) exists, and let g^ he 
defined as in [3^) . Then, for any 5 > 0, there exists a number q G M+ such that 

KAp) < -^inf (P + efi'c) < KAp) + 5 + eq, (33) 

for all r,e > 0. 

The proofs of Lemma H] and [6] are given here below. The proof of Lemma [5] follows 
from the results presented in Appendix \M This Lemma implies that p and p^ can be made 
arbitrarily close if we take e small enough. Moreover, if we express p — p'' in the anti-normal 
form, i.e., p — p^ = ^g t P'^ t^^ i^'^)* y ^^^ corresponding natural norm liip — p'") = X^st l^itl ^^ 
trivially bounded by e-e^^'^. This implies that any computer implementation of program fll9j) 
where the polynomial to minimize is expressed in normal or anti-normal form will require a 
lot of precision in order to distinguish p^ from p for low values of e. 

Proof of Lemma ^ The demonstration of this Lemma will make use of two lemmas, proven 
in Appendices [B], O, respectively. 

Lemma 7. Let s G W„ be a monomial. Then, 

t ss* I -ss* G S^. (34) 

Lemma 8. Let L be a linear functional in Wn. If L(h*h) > for any h G W„, and there 
exist c, (i > 0, fc G N such that for any monomial s G W„, the relation 

|L(.)|<dcWr(H±l±l) (35) 

holds, then there exists a non normalized quantum state p (i.e., a non-negative, trace-class 
operator) such that L{h) = tr{pTi{h)), for any polynomial h G W„. 
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Let us now proceed with the proof of Lemma HI Following Lasserre et al. [21], let yVn{r) 
be the set of elements of Wn with degree less or equal than r, and consider the semidefinite 
program: 

e* = min l{L{p)\L : >V„(2r) -> M, L linear ,L{gl) < l,L{h*h) > 0,V/i G W„,(r)}. (36) 

Noting that L = is an admissible linear functional, we have that the problem has fea- 
sible points and that e* < 0. Condition L{g^) < 1, together with Lemma [71 implies that 
the diagonal entries of all feasible moment matrices M^{y) are upper bounded, and so the 
absolute values of the rest of the entries, due to positive semidefiniteness. From these two 
observations, it follows that our problem admits a solution, i.e., e* ^ — cx3 is attainable for a 
feasible choice of L. 
The dual of (ESI) is 



max,{-e:e>0,p + ec/;'GS2}. (37) 

That this problem has solutions for any p = p* follows from the fact that, for all yU G C and 
any pair of monomials s, t G VV„, 

/is*t + /i*rs+|/i|{s*s + rt} G Si (38) 

And thus, invoking Lemma [TJ 

fis*t + I2*t*s + \l2\{ts*s t + t t*tt} G S^. (39) 

By increasing the value of e, at some point we will therefore have that p + eg"^ G S^. 

Moreover, in this particular case, there is no duality gap, i.e., the solutions of both the 
primal and dual problems coincide. Again, this can be established by invoking the quantum 
state 1^^: choosing a > such that tT{Q7i{gl)} < 1/cr, it follows that L{h) = crtr(f27r(/i)) 
is a strictly feasible point of ( 136|) and, thus, the solutions of the dual and primal problems 
are the same [20] . 

This, together with the fact that g^ — g^ is a sum of squares for s > r, implies that, for 
all e > — e* > 0, the polynomial 

P + eg' {s > r) (40) 

is also a sum of squares. 

The sequence (e*)^ is, therefore, and increasing one. We will next proof that limr_s>oo e* — 
0, and so that the e > appearing in the formulation of Lemma H] can be taken arbitrarily 
small. 

Consider the sequence L* of functionals that attain the solutions e* of the problem, and 
denote by M* their corresponding moment matrices (the entries (M*)(s j) (jj^^) where either 
||s + t||i > r or \\u + v\\i > r are assumed to be completed with zeros). By Lemma [TJ we 
have that 

(M;) (,,-),(,,-) < (M;)(,+,,o),(.+f,o) < l!i±ili±^^^dl^"+*ll^ =: d{-s, t). (41) 
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Now, perform the transformation 

{M:)^s,t),in,v) -^ {N:)i-s,F),iu,v) = d{s, t}^/^{M:)^s,t),{n,v)d{u, v)^^^. (42) 

N* is thus positive semidefinite and its diagonals are upper bounded by 1 for all r; it follows 
that all the entries of the matrices A^* are in the interval [1,-1]. By the Banach-Alaoglu 
theorem, the sequence N* admits a subsequence {rj} that converges in the weak-* topology 
to a limit A^* — )■ N when i — )■ oo [25] . Undoing the previous change of coordinates, we are left 

with an infinite sized matrix M that defines a linear functional L{p) = Yls iP{s,t}^io,s),{o,t) 
on the Weyl algebra. 

This functional satisfies L{h*h) > 0, for any polynomial h. Moreover, for any sequence 

s = S1S2, with |si| = 



1-521 



Lis)\ < (Lis,s*,)L{s;s2: 



1/2 



< 



(n+ |gi| -l)!(n+ 1^2! - 1)! 
(n-l)!(n-l)! 



1/2 



-|S1| + |S2| 



-\s\ 



< 



(n-iv: 



-r 



2n + 2 + \s\ 



(43) 



By Lemma [HI this last condition implies that there exists a non normalized quantum state 
p G SiiTi) in the Schrodinger representation such that L{h) = tr(p7r(/i)), for all h. 
Now, 



lim e* = L{p) = tr(p7r(p)) > 0, (44) 

r— >oo 

where the last inequality follows from the non-negativity assumption on p. On the other 
hand, e* < Vr, and, therefore, we have that lim^-i^oo e* — 0- 

D 



Proof of Lemma [3 Given a vector m G N", we will denote by |m) the number state 
|r7ii) ® |?Ti2) ® ... ® Ifrin)- Now, if tt{p) > 0, then Ainf(p) exists and can be written as the 
limit of a sequence of the form {(f)i\'K{p)\(j)i), where {0,} are normalized quantum states. Such 
states can be, in turn, approximated with arbitrary precision by finite linear combinations 
of number states. Choose, then, a number M ^ N such that the normalized state |$) = 
T.\\ni\\^<M dm\rh) satisfies Xini{p) < ($|vr(p)|$) < Ainf(p) + 5. 

It can be verified that, for any pair of number states \rh), \rh'), with ||m||oo, ll^'lloo < M, 
and any monomial s of the annihilation and creation operators, the inequality |(m|ss*|m')| < 



5fh fh' holds. Therefore, 



Ml 



(<l>|ss*|$) < 



(|g|+M)! 
M! 



^14 



(|s|+M)! 



M! 



(45) 



Finally, choose c > 1. It follows that 

(*l„.w) < f i{t--\t\ <'}{' + My. ^f.(u MY. 



1=0 



1)1 

(n-1)!/! 



1=0 



M\l\d 



c-1 



M+l 



(46) 



12 



where in order to identify the second and third expressions we made use of Proposition [TT] 
in Appendix |X1 Equahng to q the last result, we arrive at the promised Lemma. 

D 

Note that Lemma [H] alone provides an alternative explanation for the observed conver- 
gence to the optimal solution, this time from the point of view of the dual problem (120 p . 
The reason why /x^ does not necessarily converge in theory to Xiniip) for the exact problem 
is because the values L(s) in each program grow faster than dc^^^V {}^\^) . If, however, due 
to finite numerical precision our solvers limit the magnitude of such momenta, low order 
relaxations jjf' should provide a better approximation to AjnO 

Increasing the numerical precision of our programs should thus have a two- fold effect: 
on one hand, it should allow the computer to distinguish between p and its perturbation p. 
On the other hand, it should extend the moment matrix search space to include matrices 
with entries of very different magnitude. A high precision numerical computation should 
therefore make the curves in Figure |5] collapse to the same line. 

We used the semidefinite programming solver SDPA-GMP [U |2BJ to compute SDP ap- 
proximations of E_i,Ei with a precision of 600 digits. Figures |3] and H] show the outputs 
A*^ — A^ of both problems as a function of lambdaStar, an internal parameter of SDPA-GMP 
that constrains the magnitudes of the entries of the moment matrixo [26]. Notice that, in 
agreement with the above interpretation, the difference between A*^ and A^ = A* tends to 
zero as the constraints on the moment matrix disappear (right end of Figures E] and H]). Con- 
versely, the solutions of higher order relaxations start to differ from A^ and become closer to 
the actual solution of the problem as we restrict the magnitude of the entries of the moment 
matrix (left end of Figures [3] and S]) . 

5 Conclusion 

We have identified a paradoxical behaviour in the application to bosonic systems of the SDP 
methods widely used in quantum chemistry energy calculations and quantum information. 
Namely, we have pointed out that numerical implementations of the method seems to con- 
verge despite a simple theoretical argument showing that the first SDP relaxation should 
already provide the best lower-bound to the problem at hand. This phenomenon is similar 
to an analogous behavior observed in commutative polynomial optimization fi6\ [T7] and we 
suggested that the paradox arises from rounding errors introduced in the numerical comu- 
tation. We provided a theoretical basis for this assumption by proving that for any bosonic 
hamiltonian to be minimized there exists a perturbation of it whose ground state energy can 
accurately be approximated by the SDP method. Furthermore, we showed that the effect 
disappears as soon as we increase the computer precision. 

Our results suggest that the above problem could be avoided by constraining the values of 
the diagonal elements of the moment matrices in each program. Thanks to such constraints. 



^This does not apply to high order relaxations corresponding to /c 3> 1, as a constraint of the form 
X(o^(a*^)*) < K becomes unimplementable as soon a.s K > k\. 

•^More concretely, lambdaStar is such that (lambdaStar) • I — Mk{y) > 0. 
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Figure 3: Plot of A — A as a function of the parameter lambdaStar in the case ??i = — 1. 
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Figure 4: Plot of A — A as a function of the parameter lambdaStar in the case m = 1. 
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computer implementations of the primal problem could return reliable solutions without 
having to resort to extremely high precision numerical calculations. This approach will be 
explored in a forthcoming article. 

It is worth noting that Cimpric |2Z] proposed to use Schmiidgen's positivstellensatz for 
Weyl algebras [IB] to introduce a different SDP hierarchy than the one presented here in 
order to find rigorous lower bounds on the minimum value of arbitrary Weyl polynomials. 
The application of this method, however, requires high precision SDP solvers. 

Let us conclude with a problem for the Noncommutative Real Algebraic Geometry com- 
munity. We have shown that the set of SOS polynomials is dense in the set of positive 
semidefinite elements of the Weyl algebra, i.e., we can approximate any polynomial which is 
positive semidefinite in the Schrodinger representation by a SOS. It would be interesting to 
know if this kind of results also hold in other algebras important for quantum chemistry. For 
instance, if such an 'approximation property' were also true in algebras containing coulom- 
bian elements of the type l/|a;j — Xj\, then we would be able to estimate electronic molecular 
energies without the need of introducing orbital basis sets. 
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Appendix 

A Bounding the norm of gl 



The goal of this section is to prove the next lemma, from which Lemma [5] is a direct corollary. 
Lemma 9. Let c > 2. Then, 

higl) < ^. (47) 

The proof of this lemma, we will rely on the next two propositions. 
Proposition 10. 

^>'y = E -Im^ ^(a-)*a-. (48) 

Proof. Clearly, a*^(a*^)* = ^/^'^'.ml'^*)''^™- Evaluating the mean value of the Schrodinger 
representation of both polynomials with respect to an arbitrary coherent state \a), we have 
that 

J] .,^(a*)'a- = e-l"l^ f^ ^^^H' =: 9i\aW (49) 



? 

Z,m i=0 ■' 



It follows that 



_ ^d^g{x) 



Now, it can be proven, by induction, that 
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(50) 



x=0 



oo 



d-M = fe! ,-^ V -^^±^x^ (51) 

dx"" {k-m)\ ^j\{j + m)\ ^ ' 

ii m < k, while ^^}^ = ii m > k. The statement of the proposition follows from these 
two relations. 

D 

Proposition 11. Let |j{t G N" : ||t||i = k} denote the number of elements o/N" satisfying 
\\i\\i = k. Then, 

Proof. Our aim is to compute the number of ways in which k identical balls can be contained 
in n different boxes. Clearly, any possible configuration can be represented uniquely by a 
sequence of k dots "." and n — 1 bars "|". The number of balls ni in box 1 would then 
correspond to the number of dots on the left of the first bar; the number of balls rij inside 
box j, for 2 < j < n — 1, to the number of dots between the j — 1*'' and the j*^ bars; the 
number of balls in box n, to the number of dots on the right of the n — 1^^ box. For instance, 
the configuration ni = 2, 77-2 = 0, 71,3 = 1 would be represented by "..||.". 

It is elementary that the number of permutations oi n + k — 1 elements, out of which 
n — 1 and k are indistinguishable, is equal to (^-liiki ■ 

D 

Proof of Lemma Proposition [TU] implies that 

M«'(«'T) = t ^^^.jl^ < t ;s(i^« = 2'«- (53) 

m=0 ^ ' m=0 ^ ^ 

The last expression is logarithmically superadditive, i.e., Hj^'^'^i! < 2^A;!, for all sets of 
natural numbers {fci, ^2, •••} such that ^^ ki = k. It follows that the bound given by eq. fl53|) 
also holds for /i(a*(a*)*), with i e N", ||t||i = k. We thus have that 



^-1)! . rJrJ^*^^^U^■■¥h<k}2' 



^M) <E ,i,u.::;i^:_ .. Mavr)<E 



,_, dl*lli(n+||t||i-l)! ' ' ''-f^/n + A;-l\c'= 

\t\<r k=0 ' 

y n — 1 
k=o ^ ^ fc=o ^ ^ 



where in the third inequality we have made use of Proposition [TT] 
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B Anti-normal ordered monomials 

The following appendix establishes Lemma [3 
Proposition 12. For any /c G N, 

^fc+i(^fc+i)* _ a*ak^a'')*a e T?. (55) 

Proof. Using the OCRs, we have that 

a'a''{a^ya = -ka!"-^ {a'')* a - {k + l)a\a*)'' + a''+\a*)^+\ (56) 

Now, by induction, we have that, for any < I < k, a'' {a^)* a^~'' G S^. Indeed, for / = the 
result is obvious. Suppose now that the result holds for /. Then 

a'+^(a^)*a^-'-^ = a^a{a^ya''-^-^ = ka\a*f-^a^-^-^ + a\a^ya^-\ (57) 

and the last expression belongs to S^ by hypothesis. 
It follows that 

^k+i^j.+iy _ (fci^^aky^^ _ ^^ ^ l)al'{a^y + ka'^-^ia^ya G T?. (58) 

D 
Proposition 13. Let s G Wi he an arbitrary monomial of length k. Then, 

a\a^y - ss* G S^. (59) 

Proof. We will prove the proposition by induction. Suppose, thus, that the proposition holds 
for all monomials of length smaller or equal than k, and let s be an arbitrary monomial with 
|s| = fc + 1. There are two possibilities: 

1. s = as, with |5| = k. Then we have that 

^k+i^^k+iy _ ^^* _ a{a''{a''y - Sr)a* = J] afj*a* G S^. (60) 

i 

2. s = a*s, with \s\ = k. Then we have that 



'1 



^k+i^^k+iy _ ^^* ^ ^^k+i^^k+iy _ a*a\a''ya} + {a*{a\a''y - 5S*)a}. (61) 

The first term between brackets is a SOS by Proposition [T21 the second term is a SOS 
due to the induction hypothesis. 

To complete the induction we also have to show that the proposition also holds for A; = 1. 
But this is trivial, since, in that case, aa* — ss* equals (1), for s = a (s = a*). 

n 
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Lemma [T], Let s G Wn be a monomial. Then, 

t ss* t -ss* e si (62) 

Proof. Proposition [13] already shows that the lemma holds for n = 1. Now, suppose that 
the lemma holds for n, and let s = tu E Wn+i, with t [u) being a word with the letters 
ai, ..., an, al, ..., a* (a„+i, a*+i). Let t G N" be such that ttt*\ = a\a^)*. Then, 

t ss* I -ss* = a*aj,"|i(al,"|i)*(a*)* - a^uu*{a^)* + u | tt* | u* - uU*u*. (63) 

The first two terms on the right hand side admit a SOS decomposition due to Proposition 
[T3l The two remaining terms belong to S^ because of the induction hypothesis. D 

C States in the Schrodinger representation 

In this appendix, we demonstrate the following lemma. 

Lemma [8l Let L be a linear functional in Wn- If L[h*h) > for any h G Wn and there 
exist c, d > 0, A; G N such that for any monomial s G Wn, the relation 

|L(.)|<rfcNr(H±l±l) (64) 

holds, then there exists a non normalized quantum state p (a non-negative trace class oper- 
ator) such that 

L{h) = tr(p7r(/i)), (65) 

for any polynomial h G Wn- 

Proof. Suppose that, indeed, such a functional exists and define its characteristic function 
xiO as 



,K) = fM«_^lM=.i(e.f»«), (66) 

1=0 ^■ 

where R G W^" is the vector of polynomials 






1 
:; torm, i.e., a = ^'f^i I 

Using (jMD, we have that 



and a denotes the symplectic form, i.e., a = ©J^;^ j 



,,«),,, f(M£)M), ,a«, 



l\ 
1=0 
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with ^ = -4=maxfc{|,^2fc+i + iC,2k+2\}- That the last series converges for any value of C, follows 
from the relation 



k + 1 



This allows us to write 



a;|V"'(ix = r(^— ). (69) 



Y, j^' = / e"N|x|V^'rfx, (70) 



and the last integral converges for all a G 
Now, define the operator 



dMm^^, (71) 



- - (27r)^ 

with W^ = e*^°^^ being the so called Weyl operator |28] . 

We will next prove that, for any polynomial h, tr(p7r(/i)) = L{h). 

Because tiiW^Wfj) = {2n)''5{^ + f]), it follows that tr{pW^) = xiO = H^^)- Moreover, 
from the Weyl relations 

W^W^ = e-'^'^^/'W^^^, (72) 

it is immediate that tT{pW^Wfj) = L{e~''^'''^^'^W^_^_fj) = /(^, //). We will now show that 

f{^,v)= lim r(e,r/), (73) 

r— >oo 

Where n^,n) ^ EU=o H ^'^'f ^'"'T ) ■ 

Note that / is analytic, and that /*" corresponds to a sort of truncated Taylor expansion 
of / with respect to the variables ^, f/, where only those monomials C.'^rj^ with ||s||i, ||t||i < r 
are present. Indeed, if we take any number of derivatives of C,i,Tij (no more that r of each) 
on both sides of fl72l) and then evaluate on the point ^ = f] = 0, we will arrive at two 
polynomials pi,p2 G Wn on each side. Since this is a general relation between operators, 
both polynomials must be the same modulo the canonical commutation relations. It follows 
that 



am am' am. am' 

-————r-—f{Lv) k=f,=o= L{Pi) = L{p2) = , r{^,f]) k=n=o, 

1 1a:=1 '^SJfe 1 1/=1 CVji 1 1a;=1 "^Sifc 1 [l=l ^Vji 

(74) 
for m, m! < r. Since /(^, rj) is analytic in ^, rj, the sum of the absolute values of the coefficients 
of its Taylor expansion must converge. This implies that limr^oo f{^,v) ~ f^i^^v) ~^ fo^ 
fixed ^,f], and so we get equation f l75]) . 
Analogously, one can prove that 
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j=l \j=l 1=0 ' I 

for any k. 

Now, let V = UZiPk e Wn, with pk G {^, ^}r=i. Then, 

By taking hnear combinations of the former expectation values, we thus have that tr(p7r(/i)) = 
L{h) for any h G VV„. 

It only rests to show that p is a quantum state, i.e., it is trace-class and positive semidef- 
inite. 

From the quantum Bochner-Kinchin theorem [29], we know that xiO is the characteristic 
function of a non normalized quantum state if and only if 

i) x(0 is continuous at the origin. 

ii) for any r G N, ^i, ^2, •••, 6- ^ 1^^"; and Ci, C2, ..., c^ G C the relation 

r 

E ^'^^iXi^'^ - ^lY^'"'^''^ > (77) 

k,l=0 

holds. 

From fl68p we know that xiO is not only continuous everywhere but even analytic. On 
the other hand, from (175]) we have that 



' ^ m. — ^CfT] 



k,l=0 



^m (igfc-o--R)' 



where h^ = '^^Ck Ym^o n ■ Since hm G >V„, L{hmh'^) > 0, and so the above limit is 
non-negative. 

D 
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